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When intermediate mass stars reach their last stages of evolution they show pronounced oscillations. 
This phenomenon happens when these stars reach the so-called Asymptotic Giant Branch (AGB), 
which is a region of the Hertzsprung-Russell diagram located at about the same region of effective 
temperatures but at larger luminosities than those of regular giant stars. The period of these 
oscillations depends on the mass of the star. There is growing evidence that these oscillations are 
highly correlated with mass loss and that, as the mass loss increases, the pulsations become more 
chaotic. In this paper we study a simple oscillator which accounts for the observed properties of 
this kind of stars. This oscillator was first proposed and studied in Jl) and we extend their study 
to the region of more massive and luminous stars — the region of Super- AGB stars. The oscillator 
consists of a periodic nonlinear perturbation of a linear Hamiltonian system. The formalism of 
dynamical systems theory has been used to explore the associated Poincare map for the range of 
parameters typical of those stars. We have studied and characterized the dynamical behaviour of 
the oscillator as the parameters of the model are varied, leading us to explore a sequence of local 
and global bifurcations. Among these, a tripling bifurcation is remarkable, which allows us to show 
that the Poincare map is a nontwist area preserving map. Meandering curves, hierarchical-islands 
traps and sticky orbits also show up. We discuss the implications of the stickiness phenomenon in 
the evolution and stability of the Super-AGB stars. 

PACS numbers: 97.10.Sj, 05.45. Pq, 95.10.Fh, 82.40.Bj, 05.45.Ac 



I. INTRODUCTION 

The study of pulsating stars has attracted much atten- 
tion from astronomers. Pulsational instabilities are found 
in many phases of stellar evolution, and also for a wide 
range of stellar masses (see for an excellent review). 
Moreover, pulsational instabilities provide a unique op- 
portunity to learn about the physics of stars and to derive 
useful constraints on the stellar physical mechanisms that 
would not be accessible otherwise. Within the theory of 
stellar pulsations, there are three basic characteristics of 
the motions that the associated model may include or 
not, namely the oscillations can be linear or nonlinear, 
adiabatic or not, and radial or nonradial 0. Actual pul- 
sations of real stars must certainly involve some degree 
of nonlinearity In fact, the irregular behaviour 

observed in many variable stars is, definitely, the result 
of those nonlinear effects [pj. Moreover, the fact that 
the observed pulsation amplitudes of variable stars of a 
given type do not show huge variations from star to star 
suggests the existence of a (nonlinear) limit-cycle type 
of behaviour || . However the full set of nonlinear equa- 
tions is so complicated that there are no realistic stellar 
models for which analytic solutions exist and thus the 



investigations of nonlinear pulsations rely on numerical 
analysis. Accordingly, most recent theoretical studies of 
stellar pulsations proceed either through pure numerical 
hydrodynamical codes [p|^3[ or, conversely, they adopt a 
set of simplifying assumptions in order to be able to deal 
with a extraordinarily complex problem. Most of these 
simple models of stellar pulsation are based on a one- 
zone type of model which may be visualized as a single, 
relatively thin, spherical mass shell on top of a rigid core. 
These models have helped considerably in clarifying some 
of the complicated physics involved in stellar pulsations 
and the role played by different physical mechanisms. 

AGB stars are defined as stars that develop electron 
degenerate cores made of matter which has experienced 
complete hydrogen and helium burning, but not carbon 
burning. More luminous stars with highly evolved cores 
are called Super-AGB stars and their core is made of a 
mixture of oxygen-neon |l4|,[l5]]. It is a well known re- 
sult that the observational counterparts of Super-AGB 
stars — the so-called Long Period Variables (LPVs) - 
are radial pulsators. Moreover, long term photometry of 
these stars has shown that their light curves (the varia- 
tion of the luminosity with time) usually have some ir- 
regularities which lead to a high degree of impre- 
dictability. Due to the lack of appropriate tools for an- 
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alyzing these irregular fluctuations of the luminosity or 
the stellar radius, the scientific community did not pay 
much attention, until recently, to this category of vari- 
able stars. The development of new nonlinear time-series 
analysis tools during the last decade has changed this sit- 
uation. In particular, it has been found that these tools 
have rich applications in a broad range of astrophysi- 
cal situations, which include time analysis of gamma-ray 
bursts [ pjfl , of gravitationally lensed quasars or X- 
rays within galaxy clusters |20[ ], and of variable white 
dwarfs [plj. These tools have also been used to analyze 
theoretical models of stellar pulsations. In particular, it 
has been proven that numerical hydrodynamical mod- 
els display cascades of period doublings J^-Q as well 
as tangent bifurcations [^5 26 1. To be more specific, it 
has been confirmed |^7| that the irregular pulsations of 
W Vir models are indeed chaotic and, furthermore, it has 
been shown that the physical system generating the time 
series is equivalent to a system of three ordinary differ- 
ential equations. A similar approach has been used also 
for the study of the pulsations of other types of stars 
like, for instance two RV Tau stars, R Scuti |28| and 
AC Her (2^j, and has provided significant results con- 
cerning the underlying dynamics. This result is not a 
trivial one since it is not obvious at all that stellar pul- 
sations can be fully described by such a simple system of 
ordinary differential equations. This stems from the fact 
that the classical method to physically describe stellar 
pulsations is based on the use of a hydrodynamical code 
where the partial differential equations of fluid dynamics 
are replaced by a discrete approximation consisting of N 
mass shells. Therefore, a set of 3N coupled nonlinear 
ordinary differential equations must be solved, where N 
is typically of the order of 60. Thus, simple models do 
not only help in clarifying the basic behaviour of stellar 
pulsations but, actually, they may be able to reproduce 
with a reasonable accuracy the oscillations of real stars. 

Several such simple models have been proposed and 
studied [|30"|-[37| including the one in |Q. This model was 
intended for study the linear, adiabatic and radial pul- 
sations of AGB stars. The purpose of this paper is to 
analyze in depth the linear oscillator proposed by these 
authors, and to extend their study to more massive and 
luminous stars. The paper is organized as follows. In sec- 
tion II, we summarize the basic assumptions of the model 
and we mathematically describe the oscillator, whereas 
in section III we characterize its dynamics and describe 
a sequence of bifurcations which occur as the parameters 
of the model are varied. In section IV we compare our 
results with those for the perturbed oscillator, in order 
to get a better physical insight. In section V we discuss 
our results and, finally, in section VI we draw our con- 
clusions. 



II. DESCRIPTION OF THE MODEL 

In the search for a description which embodies the es- 
sentials of stellar oscillations, we follow the simple model 
of the driven keplerian oscillator derived in fjj], and suc- 
cessfully applied to AGB stars. The model assumes that 
the compact stellar interior is decoupled from the ex- 
tended outer layers. The driving originates in the stellar 
interior and consists of a pulsation generated by pressure 
waves. Moreover, we consider only the case of sinusoidal 
driving, where the outer layers are driven by the interior 
pressure waves that pass through a transition zone char- 
acterized by a certain coefficient of transmission. The 
motion is calculated at successive states of hydrostatic 
equilibrium. No back reaction of the outer layers on the 
inner ones is considered. 

We denominate the driving oscillator "the interior" 
and the driven oscillator "the mantle" . These are sep- 
arated by a transition zone through which the pressure 
waves from the interior propagate until they hit the man- 
tle and dissipate. The driving oscillator is represented by 
variations of the interior radius, i? c , around an equilib- 
rium position, i?o, according to: R c — Rq + aRosmu) c T 
where a and uj c are the fractional amplitude and the fre- 
quency of the driving, respectively. The mantle is rep- 
resented by a single spherical shell of mass m at instan- 
taneous radius R m . In absence of any driving force, the 
equation of motion is given by: 
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where r is the time, P is the pressure inside i? m , 
M is the mass of the rigid core. As in |5J], we assume 
that the gas follows a polytropic equation of state with 
7 = 5/3. If the region where the pulsation occurs is 
assumed to be nearly isothermal (i.e., R m f=a i?*, with 
i?* the equilibrium radius of the star), the equation of 
motion near the mantle reduces to 



dt 2 
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where we have introduced the following set of nondimen- 
sional variables: 
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FIG. 1. Poincare map for Eq.(2.3) with uj — 20.1 (typical of AGB stars) and (a) a — 0.1, (b) a — 0.4, and for uj ~ 3 (typical 
of Super-AGB stars) and fc) a — 0.1, (*d) a = 0.4. In all cases, e = 1. Note the difference in the scales of the axes. 



— cos[w(i - a; + aw 2/3 sin wi)] (2.4) 



is the characteristic dynamical frequency of the system 
and Q the transmission coefficient of the transition zone 
through which the pressure waves from the interior prop- 
agate. In the limit r w 1 (small oscillations) and absorb- 
ing some terms into t as phase shifts, one gets: 



— — x + esin[w(f — x + auj 2 / 3 sina;£)] (2.3) 

at 1 



where e = Qw ' a is the total driving amplitude, and 
x = r — 1. The parameter to is a measure of the 
core/envelope ratio, and provides information on the lo- 
cation of the source of the driving in the stellar interior. 
Note that for e = 0, Eq.(2.3) transforms into the clas- 
sical equation of the linear oscillator, x = —x. All the 
interesting features of the motion are generated by the 
perturbation (the second term of the right-hand side of 
Eq.(2.3) with e ^ 0) and its interaction with the un- 
perturbed motion. The system we are dealing with is a 
periodic-time dependent Hamiltonian system. 
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The associated Poincare map is an area preserving map. 
Hence, its dynamics does not exhibit attractors or re- 
pellers. 



III. CHARACTERIZATION OF THE 
OSCILLATOR 

The study performed in jij explored a limited set of 
values of the parameter space: e = 0.5, 0.75, 1.0, and 
a = 0.1, 0.2, 0.4 for values of uj around ~ 20, which is 
characteristic of regular AGB stars (M < 8M ). Their 
main results and conclusions are that small values of u 
produce more strongly chaotic pulsations and that for 
large values of a stable orbits in some definite regions 
of the phase space are obtained. Our aim is to extend 
the previous study to Super-AGB stars, which are char- 
acterized by a different value of uj. For the physical con- 
ditions found in Super-AGB stars, 8Af© < M < 11 M 
and i?* ~ 450 R@, it turns out that uj > 3. We will also 
explore a wider range of the parameter space than that 
investigated in jlj . Contrary to [|lj , who adopted a fourth- 
order predictor-corrector scheme for the numerical inte- 
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gration of Eq.(2.3), we use a fourth order Runge-Kutta 
integrator with step-size control and dense output as de- 
scribed in Q. We have tested several other integrators, 
specifically designed for stiff problems |^] , and we have 
obtained the same results for a given set of initial condi- 
tions. Thus we conclude that our numerical integrator is 
appropriate for problem under study. 
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FIG. 2. The effect of increasing a on the Poincare map: 
(a) a — 0.2; (b) a — 0.3 and (c) a — 0.4. e = 0.5 in all cases. 



A. AGB stars vs. Super-AGB stars: the role of ui 

In order to compare the behaviour for two different val- 
ues of uj (and, thus, the differences between pulsations of 
AGB and Super-AGB stars), Figure 1 shows the Poincare 
maps obtained with u ~ 20, typical of regular AGB stars 
(panels a and b), and those obtained for uj ~ 3 typical 
of Super-AGB stars (panels c and d), for the same initial 
conditions. As u> increases the behaviour becomes more 



irregular, and the islands and the structure of panels (a) 
and (b) quickly disappear, leading to a more chaotic be- 
haviour. Consequently, we expect that the pulsations of 
Super-AGB stars will exhibit a more chaotic behaviour 
than those of AGB stars. 

The phase portrait of our model is a typical example 
of an area-preserving map. For very small values of the 
perturbation parameter, it exhibits invariant circles on 
which the motion is quasi-periodic. Three periodic island 
chains surround the fixed point of the map. For large per- 
turbations, because of the breaking-up of the invariant 
circles and splitting of the separatrices of hyperbolic pe- 
riodic points, chaotic zones are present, as well as islands 
in the stochastic sea. Until recently, it was thought that 
the existence of islands was not very important in deter- 
mining the origin and character of chaos in Hamiltonian 
systems. However, it has been recently shown that the 
boundaries of such islands are very interesting |ic|-|4^|. In 
particular, crossing an island boundary means a transi- 
tion from a regular periodic or quasiperiodic behaviour to 
an irregular one that lies in the stochastic sea. One of the 
particularities of this singular zone is illustrated in Figure 
la. One can see a central island embedded in the domain 
of chaotic motion with a boundary separating the area of 
chaos from the region of regular motion. By changing a 
control parameter of the system, several bifurcations oc- 
cur which influence the topology of the boundary zone, 
leading to the appearance and disappearance of smaller 
islands. This may result in self-similar hierarchical struc- 
tures of islands which is crucial for understanding chaotic 
transport in the stochastic sea and the general dynamics, 
since they represent high-order resonances. 

Figure lb illustrates another essential phenomenon 
typical of some Hamiltonian systems. It is linked to the 
existence of a complex phase space topology in the neigh- 
borhood of some islands: a trajectory can spend an in- 
definitely long time in the boundary layer of the island 
(i.e., a dynamical trap). The main uncertainty in charac- 
terizing this phenomenon concerns the level of stickiness 
(that is, its characteristic trapping time) which depends 
on the parameters of the system in a still unknown man- 
ner. In general, the hierarchical structure of the islands 
in the boundary zone can explain the origin of the stick- 
iness of the trajectory within this region, which is the 
reason why this behaviour is called hierarchical-islands 
trap [^U . Much attention was generally paid to under- 
standing the structure of these islands because different 
islands correspond to different physical processes respon- 
sible for their origin. The boundary layer of an island 
may include higher-order resonant islands. Another im- 
portant reason is that the topology of the islands could 
be an indicator for the vicinity of bifurcations. However, 
no general description for the birth and collapse of hier- 
archical islands exists yet. 
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FIG. 3. Poincare maps for Eq.(2.3), a — 0.3 and different values of e. The successive birth of two dimerized island chains. 
The homoclinic and heteroclinic connections appear in red color. 



B. The role of the fractional driving amplitude, a 

In order to characterize in depth the behaviour of our 
system, we have performed a thorough parametric study 
by varying a and e, while keeping lu constant at the value 
typical of Super-AGB stars. Reasonable values of the 
fractional amplitude, a, range from 0.1 to 0.4, whereas 
the total driving amplitude, e, varies between 0.1 and 
1.0. Since a is the ratio between the amplitude of the 
internal driving and the radius of the star, the upper 
limit considered in this study is 40%, which is physically 
sound. In general, a star is characterized by its com- 
pactness, oj, and the interior-mantle coupling strength 
Q = e/(uj 4 / 3 a), which should stay within the limits of 
(0% transmission) and 1 (100% transmission) in order to 
maintain its physical meaning. Therefore some specific 
combinations of (e, a) that lay outside this range for Q 
must be regarded with caution, since they have no phys- 
ical meaning. According to all these considerations, we 
have investigated the dynamics of the system associated 
to Eq.(2.3) rewritten below in the form of a perturbed 
oscillator, 



x = y 



y = —x + esin[oj(t — x + b sin ut)] 



where b = auj^ 2 ^ , maintaining a constant value of u> ~ 3, 
as discussed above. 

In a first step, we have studied the qualitative changes 
in the dynamics (i.e., bifurcations) as the parameters e 
and a are varied. Poincare maps corresponding to the 
case a £ (0.1 — 0.4) and e = 0.5, for several initial con- 
ditions, are shown in Figure 2. The maps are character- 
ized by the same geometric structure: a region, centered 
around (x,y) = (0,0), of closed orbits, surrounded by 
a region of chaotic orbits. As it can be seen in Figure 
2, as a increases, the central region of regular behaviour 
expands to the detriment of the stochastic sea. 

We have also analyzed the behaviour of the system 
for negative values of a. A quick look at the driving 
allows one to see that negative values of a correspond 
to a mere change of phase in the driving and, therefore, 
correspond also to physically meaningful cases. However, 
and for the sake of clarity and conciseness, in this paper 
we will concentrate our efforts on the positive domain of 
this parameter. Nevertheless, we would like to remark at 
this point that the phase portrait of the Poincare map 
changes drastically depending on whether a is positive, 
negative or almost zero. 
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FIG. 4. Creation by saddle- center bifurcation of a pair of three periodic orbits for a = 0.3 and increasing values of e. (a): 
e = 0.0875 — the phase portrait at the threshold of the bifurcation; (b): the invariant circle with 6 cusps for the case (a); (c) 
e — 0.095 — two intertwined orbits of period three have been created, inside each new-born loop there is a period three elliptic 
point; (d): a sequence of heteroclinic and homoclinic connections between the new born hyperbolic points; (e): e — 0.11725 — 
the threshold of global bifurcation; (f): the new topology of the separatrices. 



C. The role of the total driving amplitude, e 

Next we have focused on values of e <C 1 in order to 
observe in detail the departure of our equation from the 
harmonic oscillator as this parameter increases. We re- 
stricted the study of the Poincare map to a rectangle 
limited by initial conditions close to reasonable values of 
the radius and velocity of the mantle. The restriction 
is also due to the approximations involved in deriving 
Eq.(2.3). For e — we get an integrable Hamiltonian 
system, whose integrals of motion are the tori given by 
the condition x 2 + y 2 = C, t = [0, 2ir/ui). The Poincare 
map has the elliptic fixed point (0,0), which corresponds 



to a periodic orbit of period P = 2ir/u> of the perturbed 
Hamiltonian system. 

We have chosen the range e € (0,0.12) and a = 0.3 
and we have obtained a cascade of local and global bi- 
furcations. Note that the considered pairs of parameters 
(e, a) are located above the curve a*(e), which is the 
curve of triplication of the elliptic fixed point. Accord- 
ing to |^3|,Q the corresponding Poincare map is in these 
cases a nontwist map. The sequence of detected local and 
global bifurcations is typical for such a map fl"5|-[48| . Fig- 
ure 3 illustrates the birth in stages of two dimerized island 
chains containing periodic points of period three. This 
is a typical scenario of creation of new orbits in nontwist 
maps. The upper-left panel exhibits the phase portrait 
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of the near-to-integrable map. The elliptic fixed point 
is surrounded by invariant circles. At e = 0.03 starts 
the birth of the first dimerized island chain, namely a 
regular invariant circle bifurcated to an invariant circle 
with three cusps. The newly born dimerized island chain 
is clearly seen in the bottom-left panel of the Figure 3, 
where a new circle with three cusps, closer to the the 
fixed point can be seen. Both dimerized island chains 
are illustrated in the bottom-right panel. Between them 
the invariant circles are meanders. Note that each cusp 
is a point of saddle-center creation |l9| . This scenario is 
a good example to illustrate the theoretical results con- 
cerning the creation of a twistless circle after triplication 

nil- 
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FIG. 5. Poincare sections for a = 0.3 and different values 
of e: (a) e = 0.2; (b) e = 0.3, and (c) e = 0.4. In red and green 
colors are represented evolutions of initial conditions close to 
the first and respectively second separatrix. 

Around the last born three-periodic dimerized island 
chain, a sequence of local and global bifurcations occurs. 



This is illustrated in Figure 4, where it can be seen that 
two independent orbits of the same period three are cre- 
ated by cusp bifurcation. They evolve in such a way that, 
finally, interact with the orbits which belong to the first 
dimerized island chain. As e increases the newly born el- 
liptic points approach the hyperbolic points of the dimer- 
ized chain. When e reaches a value of 0.11725 a global bi- 
furcation occurs: the newly created orbits interfere, and 
the hyperbolic points of the dimerized island chain be- 
come hyperbolic points with homoclinic eight-like orbits 
encircling the new created elliptic points. We have also 
noticed that as e increases, the process of creation of three 
periodic orbits repeats outside the region of the Poincare 
map we focus on, and the new orbits form chains of vor- 
tices, not dimerized island chains. This behaviour is due 
to the oscillating character of the nontwist property of 
the perturbation. 

In a second step we have paid attention to larger 
values of e. In particular we have studied the range 
e G (0.2,0.5). The increase in strength of the external 
perturbation destroys the separatrices (Figure 5, top) by 
clothing every one of them in a stochastic layer (Figure 
5, middle). As the thickness of the layers increases with 
the perturbation, depending on the positions of the sep- 
aratrices in the phase space, they can merge forming the 
stochastic sea (Figure 5, bottom). 

There are many islands, which the chaotic trajectory 
cannot penetrate. Within an island there are quasiperi- 
odic motions (invariant tori) and regions of trapped 
chaos. The stronger the chaos is, the smaller the islands 
and the larger the fraction of phase space occupied by 
the stochastic sea. The coexistence of regions of regu- 
lar dynamics (closed orbits) and regions of chaos in the 
phase space is a wonderful example of the property which 
differentiates chaotic systems from ordinary random pro- 
cesses, where no stability islands are present. This prop- 
erty makes possible the analysis of the onset of chaos and 
the appearance of minimal regions of chaos. To summa- 
rize, the system undergoes the following bifurcations, as 
it results from Figures 3, 4 and 5: 

(a) For e £ [0, 0.03), the phase portrait of the Poincare 
map is similar to that of the harmonic oscillator, 
with the elliptic fixed point (x,y) = (0,0) sur- 
rounded by almost circular closed orbits. 

(b) For e £ [0.03, 0.04), in addition to the central ellip- 
tic point, a dimerized island chain is born, namely, 
three periodic elliptic and hyperbolic points, each 
elliptic point being surrounded by a homoclinic or- 
bit to the corresponding hyperbolic point, and dis- 
tinct hyperbolic points being connected by hetero- 
clinic orbits. 

(c) For e £ [0.04, 0.07), a second dimerized island chain 
is created, which is a typical configuration after the 
triplication of an elliptic fixed point. 

(d) For e £ [0.07,0.2), a sequence of local and global 
bifurcations occurs which interphere with the first 
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(e) 



(f) 



created dimerized island chain: a regular meander 
bifurcates to a meander with 6 cusps. A new dimer- 
ized island chain is created by saddle-center bifur- 
cation. If we denote the cusps as 1, 2, . . . , 6 in clock- 
wise order, then the elliptic and hyperbolic points 
born from the odd cusps, form independent orbits 
of those born from the even cusps. This particu- 
lar dimerized island chain approaches the original 
one, and as a result the hyperbolic points of the 
two chains are connected giving rise to a pattern of 
a complicated topological type (Figure 4f). 
For e s [0.2, 0.3), we notice a more definite chaotic 
behaviour created by the instability of the remain- 
ings of the first heteroclinic orbit engulfing the pe- 
riodic points of case (c). 

For e s [0.3,0.5], we notice the extension of the 
chaotic orbits in the central region, together with 
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FIG. 6. (a): A meander near an invariant circle with 35 
cusps for e = 0.4 and a = 0.3; (b): two orbits, one filling 
the meander and the second, a chaotic orbit, confined by the 
meander. 

Another important feature of the phase portrait char- 
acteristic of nontwist maps is the existence of mean- 
ders (i.e., invariant circles which fold exactly as a me- 
ander). Meanders are created between two successively 
born dimerized island chains, and between two chains of 
vortices. In twist standard-like maps they become usual 



circles after the reconnection of the two chains. Note 
that in our case the reconnection of the first two cre- 
ated dimerized island chains, does not occur, and as e 
increases the meanders break-up leaving instead a Can- 
tor set. The top panel of Figure 6 illustrates a meander 
near a pair of dimerized island chains containing peri- 
odic orbits of period 35. Meanders appear to be robust 
invariant circles, even when the nearby orbits are chaotic 
(Figure 6, bottom panel). This behaviour was observed 
in nontwist standard- like maps p^ ] , but until now there 
is no explanation for this robustness. 



IV. COMPARISON WITH THE PERTURBED 
OSCILLATOR 

In order to get a better insight of the underlying 
physics of the oscillator studied so far, in this section we 
are going to compare it with the motion of a perturbed 
oscillator, which has been already studied extensively. 
This is important since, as it will be shown below, the 
formal appeareance of Eq.(2.3) does not differ very much 
from that of a perturbed oscillator. In order to make this 
clear consider the motion of a perturbed linear oscillator 
in the form of: 



- luqX = e sin[wx — uit] 



(4.1) 



The left-hand side of Eq.(4.1) describes very simple 
dynamics of linear oscillations of frequency loq. To ease 
the comparison with the perturbed oscillator, Eq.(2.3) 
can be rewritten as: 



where 



X + LOqX 



w(t) 



-e ra - 



•«(*)*], 



bsuiLot 
t 



1 



(4.2) 



(4.3) 



Both equations share common features. It is impor- 
tant to point out here that for our system loq = 1, as it 
results from Eq.(2.3). The Hamiltonian associated with 
the system given by Eq.(4.1) is: 



H 



(x 2 + ulx 2 ) 



■ cos[uox — Uj(t)t] 



(4.4) 



Consequently, for weak perturbations the results ob- 
tained so far should be similar to the case of a perturbed 
oscillator. And indeed this is the case, the only difference 
resides in the period 3 homoclinic orbit, that is in point 
d) of the summary given in section III.C. The creation 
of this period-3 homoclinic structure was illustrated in 
Figure 2 as a increases from a = 0.2 to a = 0.3. This 
suggests that the ultimate reason of this difference might 
be the presence of the function oj(t) from Eq.(4.3) as a 
weak detuning of resonance. 
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We rewrite Eq.(4.2) with a ^ in the following way: 

x + u x = esin u)(t — x + &sin ut)] (4-5) 

For a = 0, Eq.(4.5) transforms into: 

x + uJqX = esin[a;(t — x)] (4-6) 

All these ordinary differential equations, together with 
the well known equation for the perturbed oscillator, have 
the same mathematical structure. It is an example of 
weak chaos, where the perturbation itself creates the sep- 
aratrix network at a certain e an d then destroys it as e 
increases beyond eo by producing channels of chaotic dy- 
namics. For an unperturbed equation, the unperturbed 
Hamiltonian Hq intrinsically has separatrix structures 
and the perturbation clothes them in thin stochastic lay- 
ers, an example of strong chaos. In the phase space there 
appear invariant curves (deformed tori) embracing the 
center which do not allow diffusion in the radial direc- 
tion. Inside these cells of the web, motion occurs along 
closed-orbits, around the elliptic points from the centers 
of the cells. With the increase of e and creation of the 
stochastic layers, particles can wander along the chan- 
nels of the newly born web, a phenomenon that repre- 
sents an universal instability and gives birth to chaotic 
fluctuations. The heteroclinic structures formed as the 
perturbation increases through this bifurcation are what 
differentiates our equation from the typical equations of 
Hamiltonian chaos and in particular from those of the 
perturbed oscillator. 

V. ASTROPHYSICAL INTERPRETATION OF 
THE RESULTS 

In order to provide the reader with a feeling of the 
physical ranges for the stellar fluctuations in radius and 



velocity, in Figure 7 we show a particular time series of 
the model presented above. In this figure the variations 
of radius and velocity as functions of time are represented 
both in physical units (solar radii, km/s and years, re- 
spectively) and in adimcnsionalizcd units — as we have 
been doing so far. 

Generally speaking, the primary information that one 
can derive from a generic computed or observed time 
series is the spectral distribution of energies (or ampli- 
tudes) of the light curve. We found that for some initial 
conditions, the resultant time series do not show a single 
frequency or a small set of frequencies but, in addition, 
linear combinations of the primary frequencies may also 
appear in the spectra. This behaviour has been found 
in real stars quite frequently, and in particular in LPVs, 
Cepheids and RR Lyrae stars |50| . 

In any calculation of a time-varying phenomenon, like 
the one we are describing here, there is always the im- 
portant practical question of when the simulations should 
be terminated. For full nonlinear hydrodynamic simula- 
tions, even in the case in which a stable limit cycle seems 
to have set in, one often is worried about the fact that 
thermal changes, which take place on a longer timescales, 
could still be occurring. In a recent paper j52j, the au- 
thors have carried the calculations of Mira variable pulsa- 
tions much further than in previous investigations and, to 
their surprise, they have found that the actual behaviour 
was significantly different from what it was previously 
thought. To be precise, they obtained a new modified 
"true" limit cycle which is quite different from the ear- 
lier, "false" limit cycle. In agreement with these results, 
the numerical integration of our system presents a sim- 
ilar behavior, despite the very crude approach adopted 
here which results in an extreme simplicity of the model. 
This is shown in Figure 7 for a particular orbit: during 
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the first ~ 1400 yr the star oscillates quite regularly, and 
with small amplitudes. After this period of time, the 
quiet phase suddenly stops and pulsations become more 
violent and chaotic. Moreover, at late times the velocity 
of the outer layers exceeds the escape velocity and, hence, 
mass loss is very likely to occur, in accordance with the 
observations of LPVs, which are the observational coun- 
terparts of Super- AGB stars. 
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FIG. 8. Poincare map for the orbit of Figure 7. Panel 
(a) includes the points R(t) and V(r) for the time interval 
(1000-1400) years, whereas panel (b) shows the points R(t) 
and V(t) for the time interval (1000-2000) years. 

This behaviour — a quiet phase, followed by an ex- 
tremely violent phase — is a typical case of a dynamical 
trap discussed in section III. The Poincare map associ- 
ated with the radius and velocity variations in Figure 7 is 
shown in Figure 8. None of these figures include the os- 
cillations previous to 1000 years because the behaviour is 
similar to the one of the time interval between 1000 and 
1400 years. This example corroborates the suggestion, 
first proposed by |p3[ and later found in the very detailed 



numerical calculations of j52| that the long-term effects 
in Mira variables have an important role in our under- 
standing of the mechanism which drives mass-loss. This 
result also argues in favour of the capability of simple 
models to capture the underlying dynamics of the phys- 
ical system. To put it in a different way, our model, de- 
spite its simplicity, reproduces reasonably well the results 
obtained with full, sophisticated and nonlinear hydrody- 
namic models. It is nevertheless important to realize here 
that our model does not incorporate the secular effects 
induced by the thermal changes and, hence, it is quite 
likely that this kind of behaviour is intrinsically associ- 
ated with the physical characteristics of the oscillations 
of real stars. 

In order to get a more precise physical insight a time- 
frequency analysis would be very valuable. Consequently 
we devote the remaining of this section to such a purpose. 
Most often, astronomers have used wavelet transforms in 
their investigations of the time-frequency characteristics 
of variable light curves. However, |54| compared the re- 
sults obtained with different time-frequency methods us- 
ing both real light-curves and synthetic signals. From 
these tests they firstly concluded that the Gabor trans- 
form |5jJ provides much more informative results on the 
high frequency part of the data than the wavelet trans- 
form, but, on the contrary, they also found that the time- 
frequency analysis using the Choi- Williams distribution 
p6[ is definitely superior to both methods (at least on 
the data they used). Further investigations carried out 
by the same authors |HJ led them to conclude that in 
general one cannot claim a priori that any of these meth- 
ods is better than the others. In fact, this depends largely 
on the nature of the signal and on which kind of features 
one tries to enhance. Therefore, it is always advanta- 
geous to use simultanously several of them. However, as 
the Choi- Williams distribution is generally accepted to 
be the most useful among all the methods, we will center 
our analysis using this tool, as it was done in |plj| . 

The time- frequency analysis of a model computed with 
lu ~ 3, a — 0.3, e = 0.5 and the initial condition 
(xo,yo) — (0.0,0.02), is shown in the central panel of 
Figure 9. In order to be used as a visual help this time 
series has been reproduced in the bottom panel of this 
figure. We have split the time series into two pieces, in 
order to better illustrate the stickiness of the oscillator, 
and we have computed the Fourier transforms of both 
pieces separately. These Fourier transforms are shown 
in the top panels of Figure 9, the left one corresponds 
to t < 1400 yr, whereas the right one corresponds to 
t > 1400 yr. Note the difference in the scales of the 
power, which is much larger for the right panel. As it 
can be seen in these panels the power is suddenly shifted 
from relatively large frequencies to much smaller frequen- 
cies. 
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and for r > 1440 yr (top tight panel), and the corresponding time series (bottom panel). 



This behavior is even more clearly shown in the Choi- 
Williams distribution. Note as well the contribution of 
the frequency / ~ 2.1 yr -1 at small times, just before 
the beginning of the burst (at r < 1430 yr). Indeed, this 
small contribution could be at the origin of the transfer 
of power to / ~ 1.2 yr" 1 , which occurs inmediately after 
t ~ 1430 yr and which ultimately leads to the burst at 
t ~ 1440 yr. Also remarkable is the small time elapsed 
since the beginning of the transfer of power, which is only 
about 10 yr. 

VI. CONCLUSIONS 

We have presented and analyzed the complex dynam- 
ics of a forced oscillator which is interesting not only from 
the mathematical point of view, but also because it de- 
scribes with a reasonable degree of accuracy the main 
characteristics of some stellar oscillations. This model 



has been previously used to study the irregular pulsa- 
tions of Asymptotic Giant Branch stars |jj , and we have 
used it to study the pulsations of more massive and lumi- 
nous stars, the so-called Super- Asymptotic Giant Branch 
stars. In doing this we have extended the previous stud- 
ies to a range of the parameters specific for this stel- 
lar evolutionary phase. We have found, in agreement 
with the observations of Long Period Variables which 
are the observational counterparts of Super- Asymptotic 
Giant Branch stars, that the oscillator shows a chaotic 
behavior. It is important to realize as well that this kind 
of stars shows a more pronounced chaotic behavior than 
regular Asymptotic Giant Branch stars of smaller mass 
and luminosity. We have also characterized in depth the 
full sequence of bifurcations as the physical parameters 
of the model are varied. We have found a rich set of 
local and global bifurcations which were not described 
in [Q. Among these, perhaps the most important one 
is a tripling bifurcation, but meandering curves, hierar- 
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chic al islands traps and sticky orbits also appear. Cor- 



ings of the 6 th Experimental Chaos Conference, astro- 



res ppndingly, the resulting time scries also show a rich 
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behaviour. In particular, we have found that although 
there are light curves which show a rather regular be- 
haviour for certain values of the parameters of the phys- 
ical system and given initial conditions, there are as well 
some other light curves which show clear beatings or lin- 
ear combinations of two main frequencies up to terms of 
2/o + 7/i, being / the fundamental frequency and fx 
that of the first overtone, and even more complex or- 
bits. For the parameters and initial conditions leading 
to more irregular behaviour, we noticed the existence of 
both clear chaotic pulsations and sudden changes from 
a limit-cycle to chaotic pulsations, the latter being asso- 
ciated with the stickiness phenomenon characteristic of 
some Hamiltonian systems. For these orbits the velocity 
of the very outer layers clearly exceeds the escape veloc- 
ity. Hence, for these chaotic pulsations mass-loss is very 
likely to occur, in good agreement with the observations, 
which correlate the degree of irregularity with the mass- 
loss rate. Regarding the stickiness of some orbits which 
we have found for a set of parameters, perhaps the most 
important result is that the long-term effects found in real 
stars are reproduced by our model, despite of its simplic- 
ity, even though the driven oscillator studied here does 
not incorporate the effects of secular changes. Hence it 
is quite likely that this kind of behaviour which has been 
already found in full hydrodynamical simulations [ p2"| , is 
intrinsically associated with the physical characteristics 
of the oscillations of real stars and not to the long-term 
thermal changes. 
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